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ABSTRACT 

We use several quasar samples (LBQS, HBQS, Durham/AAT and EQS) to determine 
the density and luminosity evolution of quasars. Combining these different samples and 
accounting for varying selection criteria require tests of correlation and the determina- 
tion of density functions for multiply truncated data. We describe new non-parametric 
techniques for accomplishing these tasks, which have been developed recently by Efron 
and Petrosian (1998). With these methods, the luminosity evolution can be found 
through an investigation of the correlation of the bivariate distribution of luminosities 
and redshifts. We use matter dominated cosmological models with either zero cosmo- 
logical constant or zero spatial curvature to determine luminosities from fluxes. Of the 
two most commonly used models for luminosity evolution, L oc e kt ^ and L oc (1 + z) k ' , 
we find that the second form is a better description of the data at all luminosities; we 
find k' = 2.58 ([2.14, 2.91] one a region) for the Einstein - de Sitter cosmological model. 

Using this form of luminosity evolution we determine a global luminosity function 
and the evolution of the co-moving density for the two types of cosmological models. 
For the Einstein - de Sitter cosmological model we find a relatively strong increase in 
co- moving density up to z ^ 2, at which point the density peaks and begins to decrease 
rapidly. This is in agreement with results from high redshift surveys. We find some co- 
moving density evolution for all cosmological models, with the rate of evolution lower 
for models with lower matter density. We find that the local cumulative luminosity 
function <&(L ) exhibits the usual double power law behavior. The luminosity density 
C(z) = / °° L^(L, z)dL, where ^(L,z) is the luminosity function, is found to increase 
rapidly at low redshift and to reach a peak at around z ~ 2. Our results for C(z) are 
compared to results from high redshift surveys and to the variation of the star formation 
rate with redshift. 
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1. INTRODUCTION 



Investigations of the evolution of the quasar population have played a major role in the devel- 
opment of our ideas about the nature of these sources and their connection to other extragalactic 
objects. Ever since the first complete survey of 3C radio quasars by Schmidt (1968) and the sub- 
sequent survey of 4C quasars by Lynds and Wills (1972) it has been evident that the population 
of quasars as a whole has undergone rapid evolution. Using the so-called (V/V max ) method these 
authors interpreted the evolution with redshift z as caused by an increase in the co-moving density 
of quasars with redshift. However, both the source counts (Giaconni et al. 1979, Tananbaum et al. 
1979) and the redshift distribution of optically selected samples of quasars (see e.g. Marshall 1985) 
clearly showed that such pure density evolution (PDE) models, for which the luminosity function 
is separable as 

V(L,z)=iP(L)p(z), (1) 

cannot be correct. As more data was accumulated the pure luminosity evolution (PLE) model, 
with 

*(L,z) = i/>(L/g(z))/g(z), (2) 

gained more popularity. The function g(z) describes the luminosity evolution of the population 
and L Q = L/g(z) is the luminosity adjusted to its present epoch (z = 0) value. This model, while 
providing a better fit to the data than that of pure density evolution, also appears to be inadequate 
in many cases (see e.g. Petrosian 1973, Schmidt and Green 1978, Koo and Kron 1988, Caditz and 
Petrosian 1990). Without loss of generality, we can write the luminosity function as 

*(L,z) = p(z)il>(L/g(z),a i )/g(z). (3) 

With ip normalized such that f£° ip(L,ai)dL = 1, p(z) gives the density and its evolution and 
vp(L ,ai) describes the local luminosity function (with g(0) = 1). Here we explicitly include the 
«i, parameters such as spectral index and break luminosity that describe the shape of the lumi- 
nosity function. In general, these parameters may vary with redshift. A surprising, and a priori 
unexpected, result has been the absence of evidence for strong shape variation. 

Such results imply that the density and luminosity functions p(z) and g(z) describe the physical 
evolution of sources; e.g. the rate of birth and death of sources and the changes in source luminosity 
with time. Cavaliere and colleagues (see Cavaliere and Padovani 1988 and references cited therein) 
were the first to emphasize this fact. A more complete description of the relation between the 
physical evolution and the functions describing the generalized luminosity function (or statistical 
evolution) can be found in Caditz and Petrosian (1990). Unfortunately, this relation is not unique, 
thus any test of quasar models via their expected evolution with cosmic epoch must usually involve 
additional assumptions. For example, quasars could be long lived (compared to the Hubble time) 
sources created during a relatively short period at high redshift undergoing continuous luminosity 
evolution. This is the model used in Caditz, Petrosian and Wandel (1991). Alternatively, quasars 
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could be short lived phenomena with birth rate, death rate and luminosity that vary systematically 
with redshift (see Siemiginowska and Elvis 1997). Most of the work along these lines assumes what 
is now the standard model for the source of energy of quasars and other active galactic nuclei (or 
AGNs); namely, accretion onto a massive black hole. 

Another interesting aspect of quasar evolution is the relation between the evolution of the 
luminosity density of quasars, £ = / L^(L, z)dL, and similar functions describing the evolution of 
galaxies, such as the star formation rate (SFR). As shown by the high redshift surveys of Schmidt 
et al. (1995) and others (e.g. Warren et al. 1994), the rapid rise with redshift of the luminosity 
of quasars stops around redshift 2 or 3 and is expected to drop, perhaps mimicking the SFR (see 
Shaver et al. 1998, Cavaliere and Vitorini 1998, and Hawkins and Veron 1996). 

In this paper we determine the luminosity function and its evolution for combined samples 
of quasars from various surveys, described in §2. For a complete review of the various ways of 
accomplishing this task, see Petrosian (1992). We use here non-parametric methods based on 
Lynden-Bell's (1971) idea which was generalized by Efron and Petrosian (1992) to allow not only 
a determination of the functional forms but also a test of correlation between redshifts and lumi- 
nosities. This is essential for any determination of the functional form of the luminosity evolution 
g(z). The above papers deal with magnitude limited samples with only an upper magnitude (lower 
flux) limit. However, because of observational constraints, some of the samples are truncated in 
redshift space and some have an upper magnitude (lower flux) limit. The methods described in 
the above papers cannot be used for such multiply truncated data. New methods for treating this 
type of data were developed by Efron and Petrosian (1998). In §3 we describe these new techniques 
and their relation to the older methods. The choice of cosmological model also plays a crucial role 
in such determinations. It is well known that one can not determine both the evolution of the 
luminosity function and the parameters of the cosmological model from magnitude limited samples 
alone. Only by assuming values for the cosmological parameters such as matter density, curvature 
and cosmological constant can one calculate the form of *&(L, z). Alternatively, with some assump- 
tions about one can test various cosmological models. In §4 we describe the cosmological 
model parameters and in §5 we present the results of applications of the new statistical methods 
to the data described in §2. Finally, in §6 we summarize these results and present our conclusions. 

2. THE DATA 

There have been many surveys of quasars, ranging from the original radio selected samples of 
the 3C (Schmidt 1963) and 4C (Lynds and Wills 1972) surveys to a variety of other optically and 
X-ray selected samples. In this paper, we will focus only on optically selected data. We use data 
from four samples with relatively similar selection criteria that provide a somewhat homogeneous 
data set spanning a large area of the L — z plane. In order to combine the samples, all magnitudes 
were transformed to the B band. Some of the samples were artificially truncated in order to insure 
completeness within the truncation limits. The combined sample consisted of 1552 objects in the 
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range 15.5 <J B <; 21.2 and 0.3 < z < 2.2, with well defined upper and lower magnitude limits for 
each object. 

Figure 1 shows the distribution of the complete surveys in the B — z plane, which at first sight 
shows little evidence for a Hubble relation of B oc 51og(z). However, as we shall discuss below, 
there is evidence for cosmological dimming with redshift of the sources. 
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Fig. 1. — The B magnitude - redshift data for the four samples LBQS, HBQS, Durham/ AAT and 
EQS described in §2, for 0.3 < z < 2.2. The EQS data are marked with plus signs. The data 
are fitted to the parametric form B{z) = B — (3\og(d 2 L (z)K(z))+ constant, where di{z) is the 
luminosity distance and K(z) is the K-correction term. The best fit {fi = 0.84 for the Einstein - de 
Sitter cosmological model) is shown by the dashed line. The solid line shows the expected relation 
between B and z for standard candle sources without luminosity evolution ((3 = 2.5) in the Einstein 
- de Sitter model. 

2.1. The Large Bright QSO Survey 

The Large Bright QSO survey (LBQS) contains 1055 QSOs in the magnitude range 16.0 < 
Bj < 18.85 and redshift range 0.2 < z < 3.4 (Hewett et al. 1995). The faint magnitude limits for 
the 18 fields range from 18.41 to 18.85 and the bright magnitude limit is 16.0 for the entire sample. 
In order to transform from Bj to B magnitudes the color equation of Blair and Gilmore (1982) 




B = Bj + 0.28(5 - V) 



(4) 
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was used assuming an average B — V of 0.3 for the entire sample. When we combine this sample 
with the AAT sample and others we introduce artificial cutoffs at z = 0.3 and z = 2.2, which 
are the completeness limits for the AAT data. In addition, we removed all objects brighter than 
B = 16.5 in order to insure completeness near the bright magnitude limit, leaving a total of 871 
objects. 

2.2. The Homogeneous Bright QSO Survey 

Data from the six deepest fields of the Homogeneous Bright QSO Survey (HBQS) were pub- 
lished by Cristiani et al. (1995) The six deepest fields of the HBQS contain either Bj or B' 
magnitudes for 285 QSOs in the range 15.5 < Bj < 18.85, with faint magnitude limits ranging 
from 18.25 to 18.85. The bright magnitude limits are generally lower than those of the LBQS, vary- 
ing from 14.0 to 16.0. For observations in the Bj band equation (||) was used and for observations 
in the B' band the equation 

B = B' + 0.11(5 -V) (5) 

of Blair and Gilmore (1982) was used assuming B — V = 0.3. Again, when combining this sample 
only objects in the redshift range 0.3 < z < 2.2 were used, leaving a total of 254 objects. 

2.3. The Durham/ A AT Survey 

The Durham/AAT survey contains 419 QSOs in the magnitude range 17 < b < 21.27 (Boyle et 
al. 1990) and redshift range 0.3 < z < 2.2, giving information about the QSO luminosity function 
in a different regime than the previous two samples. The faint magnitude limit ranges from 20.25 
to 21.27 and the bright magnitude limit ranges from 16.4 to 18.0. The data are given in the "6" 
magnitude system, which according to Boyle et al. (1990) may be converted into the B system by 
the relation 

5 = 6 + 0.23(5-^-0.9). (6) 
As above, the average value B — V = 0.3 was used to determine B. 

2.4. The Edinburgh QSO Survey 

A subsample of the Edinburgh QSO Survey (EQS) consisting of 12 QSOs brighter than B = 
16.5 was published by Goldschmidt et al. in 1992. Of these, the 8 that fall in the redshift range 
0.3 < z < 2.2 were added to the combined samples to give information about the luminosity 
function at the bright end. 

There have been several previous analyses of the above quasar surveys. Boyle et al. (1990) 
used binning techniques to fit the the Durham/AAT data to the PLE model. More recently, La 
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Franca and Cristiani (1996) fit the LBQS and HBQS data to a more complex luminosity function 
(involving 3 or more parameters) with no density evolution. Hatziminaoglou et al. (1998) examined 
the PLE and PDE cases separately using both the Durham/AAT and LBQS data. As described 
below, we combine all of these data and determine both the luminosity evolution g(z) and the 
density evolution p(z). 

3. STATISTICAL METHODS 

The statistical problem at hand is the determination of the true distribution of luminosities 
and redshifts of the sources from a biased or truncated data set, such as a flux limited sample. 
Several different techniques exist that give such a determination. The most common method is to 
bin the data and fit to some parametric forms of ip(L) or p(z). However, it is preferable to avoid 
binning and to use non-parametric methods whenever possible. For a review of the various methods 
see Petrosian (1992). Assuming the general form of equation (||) for the luminosity function we 
must determine the functional forms of the local luminosity function -0(L o ,aj), density evolution 
p(z) and the luminosity evolution g{z) as well as the changes in the parameters on with redshift, 
if any. This last aspect is a higher order effect and will not be within the scope of this paper. We 
will discuss the certainty with which this can be ignored in the final analysis. 

All non-parametric techniques for determining the distribution in a bivariate setting require 
that the data be expressed in terms of two uncorrelated variables, i.e. that we use variables x 
and y for which the density function is separable: ^(x, y) = p{x)%l){y). Thus before applying non- 
parametric methods one must first determine the degree of correlation of the data in the x — y plane. 
This determination and the process of removing the correlation is equivalent to a determination 
of the functional form of the luminosity evolution g{z) and the subsequent transformation L — > 
L Q = L/g{z). This cannot be accomplished easily by non-parametric methods. Therefore, we chose 
two parametric forms for the luminosity evolution, gk(z) and gk'(z), and find the values of the 
parameters k and k' for which L Q and z are uncorrelated. Once this is done, the non-parametric 
methods described in Petrosian (1992) may be used to determine p{z) and tp(L ). 

The methods normally used to test correlation and determine the distributions are suited 
for simple truncations, such as y < f(x) (which by defining x' = f(x) can be reduced to the 
generic case of y < x'). This is sufficient for simple flux limited data. However, the majority of 
astronomical data, and quasar data in particular, suffer from more than one truncation. The data 
may have an upper as well as a lower truncation, f~(x) < y < / + (x), or there may be similar 
truncations in the value of the other variable. In addition, the functions f~(x) and f + (x) may 
not be continuous or even single valued. In general, multiply truncated redshift - magnitude data 
may be written as {zj, m«, [z~ , zf], [m~ , mf\]f =l where [z£ , zf] and [m^ ,mf] are the observational 
limits on z and m for the i th object, respectively. Given a cosmological model £1 this gives data 
of the form [zf,^], [L^ , Lf\\f =l . If one assumes that each object has the same redshift 

limits [z~ ,zf], which is the case in the majority of our analyses, then the problem is to test the 
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correlation and distribution of x and y from a data set {xi,yi}^L 1 given truncation limits [y^ ,yf\ 
for each point. The previous methods developed for this test (see Petrosian 1992 and Efron and 
Petrosian 1992) are suited for one sided truncations. In a more recent work (Efron and Petrosian 
1998) we have developed methods, which are a generalization of the earlier methods, for dealing 
with doubly truncated data. We will briefly review these new methods of testing correlation and 
non-parametrically determining the density evolution and luminosity function. 



3.1. Tests of Correlation and Determination of Luminosity Evolution 



3.1.1. Untruncated data 



If x and y are independent then the rank Ri of Xi in an untruncated sample (i.e. a sample 
truncated parallel to the x and y axes, so that yf are independent of Xj) will be distributed uniformly 
between 1 and N with an expected mean E = \(N + 1) and variance V = ±(N 2 - 1). One may 
then normalize Ri to have a mean of and a variance of 1 by defining the statistic Tj = (Ri — E) /V . 
One then rejects or accepts the hypothesis of independence based on the distribution of the Tj. 

One simple way of doing so is by defining a single statistic tdata based on the T« with a mean 
of and a variance of 1. One then rejects the hypothesis of independence if \tdata\ is too large (e.g. 
\tdata\ > 1 for rejection of independence at the 1 a level). The quantity 

is one choice of such a test statistic. This r is equivalent to Kendell's r statistic, which is defined 
in the following manner. Consider all possible pairings V = between data points and call a 

pairing positive if (xi — Xj)(yi — yj) > and negative if (xi — Xj)(yi — yj) < 0. If there are no 
ties then the r of equation (^) is equivalent to Kendell's r statistic 

# positive G V — # negative G V 



3.1.2. Data with One-Sided Truncation 



(8) 



A straightforward application of this method to truncated data will clearly give a false cor- 
relation signal. Efron and Petrosian (1992), and independently Tsai (1990), describe how this 
method can be applied to data with one-sided truncation, i.e. either y^ = — oo, % = 1,... ,N or 
yf = do, i = 1, . . . , N. For example, if yf = oo but yj varies with x^ then the above procedure is 
modified as follows. For each object define the comparable or associated set 



Ji = {j ■ Vj > Vi, Vj < Vi} 



(9) 
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to consist of all objects of greater y for which the value y = y% could possibly be observed. This 
is the same set as defined in Lynden-Bell's C~ method and, unlike the definition of Efron and 
Petrosian (1992, e.q. (2.9)), does not include the object in question. In the case of luminosity- 
redshift data this is the largest subset of luminosity (not magnitude or flux) and volume limited 
data that can be constructed for each given (Li, z{). If x and y are independent then we expect the 
rank Ri of Xi in the eligible set, 

Ri = #{j G Ji ■ Xj < Xi}, (10) 

to be uniformly distributed between 1 and N%, where iV, is the number of points in Jj. The rest 
of the procedure follows as for untruncated data. The normalized statistic Tj is defined here as 
Ti = (Ri - Ei)/Vi where Ei = \(Ni + 1) and Vi = ^{Nf - 1). The test statistic r is then defined 
by 

and is equivalent to a version of Kendell's r statistic defined by equation (||). In this case, however, 
we consider only the set of possible pairings allowed by truncation V = {(i,j) : y% > yJ,Uj > U^}- 



3.1.3. Multiply Truncated Data 



A generalization of the above method to doubly (or multiply) truncated data was developed 
recently by Efron and Petrosian (1998). The method is equivalent to the previous method, with 
the eligible set defined as 

Ji = {j ■ Vj > Vi^ Vi e (yj,yf)} (12) 

and the set of allowed pairings 

V = {(i,j) : y % g (yj,y+), yj e (vi ,yf)} (13) 
defined such that each object lies within the truncation region of the other. 

In this case, however, the distribution of the rankings (or of r) is unknown. If the data are 
uncorrelated then r must still have a mean of zero and a bootstrap method may be used to determine 
the variance V T as follows. By assuming the data {xi, yi} are uncorrelated one can use the methods 
of the next section to determine ip(y), the underlying (i.e. non-truncated) probability distribution 
of y. Once ip(y) is found one can simulate N s i m sets of data with underlying probability density 
fp(y) and truncation limits [y~ ,yf] for the i th object in each set. For each simulated set of data 
Dfc one may find as in equation (0) and estimate V T from the distribution of {Tfc}^™- For large 
numbers of simulations N s i m the error in this determination of V T is approximately V T /V ^sim ■ 
Given V T one can define a normalized test statistic t/V t with a mean of and a variance of 1. 
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3.1.4- The Luminosity Evolution 

If x and y, in this case the luminosity and redshift, prove to be independent, which would 
be the case if \t\ < 1, one may assume that there is no luminosity evolution and proceed with 
the determination of the univariate distributions ip(L) and p(z) of equation (]!]) using the methods 
described below. However, if |r| > 1 then L and z cannot be considered independent and one may 
assume that the most likely explanation is the presence of luminosity evolution (g(z) ^ constant). 
Another possibility is the variation of the shape parameters ctj with z. We will return to this 
possibility below. One can determine the function g(z) parametrically as follows. 

Given a parametric form for luminosity evolution Qk{z) one can transform the luminosities into 
L Q {k) = L/gk{z) and proceed with the determination of the correlation r(/c) between L Q and z as a 
function of the parameter k. If r is normalized to have a standard deviation of 1, then the values of 
k allowed by the data at the 1 a confidence level are {k : \r(k)\ < 1}, and the most likely values of 
k are those with r(k) = 0. Although in principle it is possible for the function r(k) to have several 
zeroes, this did not occur for the specific cases described in §5. 

3.2. Non-Parametric Determination of Distribution Functions 

Once a function g is found that removes the correlation between x and y (or z and L), the task 
is to find the underlying distributions p(x) (the density evolution) and ip{y') (the local luminosity 
function) given uncorrelated data {xi,y[ = y%/ g(xi)}f =1 and truncation limits If the 

truncation is one-sided (e.g. y- + = oo for all i) then a variety of non-parametric methods can 
be used to determine the univariate distribution functions. As shown by Petrosian (1992) all 
non-parametric methods reduce to Lynden-BelFs (1971) method. As with the tests of correlation 
described above, the gist of this method is to find for each point the comparable set defined in 
equation (||) and the number N{ of points in this set. For example, the cumulative distribution in 
y', $(y') = /J 3 ip(t)dt, is given by 




For doubly truncated data the comparable set is not completely observed, thus a simple analytic 
method such as the one described above is not possible. However, it turns out that a simple iterative 
procedure can lead to a maximum likelihood estimate of the distributions. Here we give a brief 
description of this method; for details the reader is referred to Efron and Petrosian (1998). 

Assume that the underlying density function ifj(y) is discretely distributed over the iV observed 
values of y. If we let tpi = tp(yi) be the probability density at yi then $j = Y^j^Pji where the 
summation includes all data points for which yj S [y^,Vi], is the total probability density for the 
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truncation region [y, ,yf\. If we define the matrix J by 

fl, if Uj £ []Ji~, yf\ , , 

v ~\o, it yj t\yr,vt] (5J 

then the definition of $j is equivalent to = J • Vl/ where SI/ = (ipi, . . . , i^n) arid = ($1, . . . , 3>iv)- 
This matrix J contains all of the information about the data and the truncation limits needed to 
find the vector of probability densities Vl/ . We also have the normalization condition on ¥, 

N 

E^ = l- (16) 

8=1 

From the definition of it follows that the conditional probability ip(yi\ [y~ , y^]) of observing a 



value yi within the truncation region [y- ,yj] is 



MvWv- v + ]) ~ if Vi G [Uj ,V ^ ] (17) 

l{!l ' !j i-vj> if//, - //,.//;;• (l7) 



The final condition on Vl/ is determined by maximizing the likelihood of observing the actual data, 

N N i 

p d ata = n^ifor = n ^ t / • (i8) 

i=i i=i 

By setting dP ^ ta = it follows that t/^ 1 = J2j Jjk^j 1 ,^ = 1,...,N. This may be written 
compactly as 

i=jt_L_ (19) 

with the notation ^ = (a± , . . . , aj^ 1 ). Thus we have reduced the problem of finding the density 
function for data with arbitrary truncation to the "moral equivalent" of an eigenvalue problem for 
the matrix J. 

In practice, this condition may be used as a recursive formula to determine Vl/. One starts with 
an initial guess for the density vector Vl/ . Equation ( |l9|) then gives the recursion relation 

' J^ETJT+c^ (20) 



where the constant is determined by the normalization condition Y^i^i = 1- One m ay use 
as an initial guess the untruncated solution ip° = In most problems, however, one of the two 
truncations will have a more pronounced effect than the other. In this case, one may ignore the 
weaker truncation and use the result based on the one-sided method as an an initial guess. We 
found that with this initial guess and data confined to one region of the L — z plane the sequence of 
M/^') defined by equation (|20j ) usually converged quickly. However, for combined samples spanning 
different regions of the L — z plane it was helpful to use an algorithm to accelerate the convergence of 
the series of For this purpose we used Aitken's 5 2 method, which gives an improved estimate 

for the terms of series by assuming approximately geometric convergence (see, e.g. Press et al. 
1992, chapter 5.1). 
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4. COSMOLOGY AND MODELS OF LUMINOSITY EVOLUTION 



In order to determine the intrinsic parameters of an object from the observed data one must 
assume a certain cosmological model. Many cosmological models may be described in terms of a few 
fundamental parameters, which (see e.g. Peebles 1993) for a matter dominated (non-relativistic) 
universe are the matter density p , the cosmological constant A and the curvature of space k (which 
is +1, or —1 for closed, flat or open universes, respectively). Using Hubble's constant H Q these 
parameters may be written in dimensionless form as 

_ 8ttGp kc 2 , A 

= w k = ~wm" n " = m (21) 

where R Q is the value of the expansion parameter of the universe at the present epoch. These 
three parameters are related via the Friedman-Lemaitre equation S]« + + Qa = 1, allowing us 
to eliminate one of them in favor of Hubble's constant. For example, the curvature term may be 
written 

n k = i - (n M + n A ). (22) 

We will consider the two classes of cosmological models given by f^A = (no cosmological constant) 
and Qfc = (flat universe with cosmological constant). For calculations of the luminosity function 
we will pay particular attention to the two cases £lk = f^A = and Qa = 0.85, SIm = 0.15. The 
first of these is the standard Einstein - de Sitter model, and the second is an inflationary model 
with parameters in accordance with current observations. For definiteness, Hubble's constant was 
assumed to be H Q = 70 km/(s Mpc), although most results are independent of this assumption. 

Once the values of the cosmological parameters are fixed, calculations of intrinsic parameters 
are relatively straightforward (Peebles 1993). The absolute luminosity, for example, takes the form 

L = f^d\K{z) (23) 

where / is the observed flux, the luminosity distance di is 

c ,„ . sin Vku(z) , 4 . 



and K(z) is the K-correction term. The co-moving coordinate distance is 

dz 

lo ^n M (i + zf + n k (i + z)' 2 + n A 

and the co-moving volume contained within a sphere of radius corresponding to redshift z is 



r z dz 

«(*) = / ; 7 , s., = (25) 

V ; Jo x/Vm l + z 3 + 17, 1 + z ) 2 + Qx 



r t z dii d 2 

Viz) = 4tt(— ) / ~ ^ dz. (26) 
V ; K H J J dz{l + z) 2 V ; 



In general, these integrals must be evaluated numerically. 
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In order to determine the K-correction term K(z) one must make an assumption about the 
quasar spectrum in the optical region. The general practice here is to assume a power law spectrum 
^optical « 1/01 with spectral index a ~ —0.5, which gives K[z) = (1 + z) 1+a ~ + z). 

Two models for luminosity evolution were used: gk{z) oc e kt<yZ ^ and gk(z) oc (1 + z) k . The first 
of these assumes an exponential dependence on the fractional lookback time t(z), which is defined 
as t(z) = 1 — T(z)/T(0), where T(z) is the age of the universe at redshift z, 

x 1 f°° du dz , . 

TM-^'jf dzJTTT)' <2T> 

The second model assumes a power law dependence on the scale factor of the universe (or the 
expansion parameter R), which is independent of the cosmological parameters. Analyses of earlier 
data (see e.g. Caditz and Petrosian 1990) have traditionally given estimates of quasar evolution 
for these two parametric forms of g{z) « e 7 - 5t ( z ) an d g(z) ~ (1 + z) 3 . 



5. THE EVOLUTION OF THE LUMINOSITY FUNCTION 



In what follows we apply the procedures of §3 to determine the correlation between the lumi- 
nosities and redshifts and test parametric forms for the evolution of the luminosity, the function 
g(z). Then we transform all luminosities to their present epoch values L Q = L/g(z) and deter- 
mine the co-moving density evolution p(z) and the present epoch luminosity function ip(L ) for the 
cosmological models described in §4. We apply these tests to the surveys described in §2 individ- 
ually and in various combinations. Before presenting these results we discuss briefly the redshift - 
magnitude data, i.e. the Hubble diagram shown in Figure 1. 



5.1. The Quasar Hubble Diagram 



As is evident from Figure 1, at first glance there seems to be very little correlation between 
the redshifts and magnitudes (or fluxes) of quasars; i.e. there is no obvious evidence for a Hubble 
type relation. This result is well known. For example, a preliminary test of correlation between m 
and z in a small subsample by Efron and Petrosian (1992, Figure 6) showed no correlation. Earlier, 
the absence of a clear Hubble relation was used as an argument against the cosmological origin of 
quasar redshift (see, e.g. Burbidge and O'Dell 1973). This is not the only possible interpretation, 
however. One would not expect a simple Hubble diagram for sources with a broad luminosity 
function (non-standard candle sources, see e.g. Petrosian 1974). The absence of an obvious Hubble 
relation can also arise from approximate cancelation between cosmological dimming and luminosity 
evolution. Exact cancelation of these two effects is highly implausible and could bring into question 
the basic assumptions about the distribution of sources. To clarify this situation we have applied the 
correlation tests described in §3 to the surveys of §2. First ignoring the high flux (low magnitude) 
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limit, i.e. with m m i n = —oo, we use the one-sided tests and find the results labeled t\ shown 
in Table 1. When using the double-sided tests we find the results T2, which indicate slightly less 
correlation, as expected because the one-sided methods ignore the slight truncation induced by the 
high flux limits. These tests, when applied to the combined data sets, give a correlation of r = 3.63. 
This result rejects the hypothesis of independence between B and z at the 99.97% confidence level 
and is independent of any cosmological parameters. In addition, we may test the parametric fit 
B(z) = B — f3\og{d\(z, Q)K(z))+ constant for the data using the methods for multiple truncations 
to determine the best value of f3, i.e. the value for which B{z) and z are uncorrelated (r = 0). 
A value of f3 = 2.5 is what one would expect for standard candle sources with a very narrow 
luminosity function, while a value of (3 = would mean the complete absence of a Hubble relation 
and the exact cancelation described above. The results shown in Table 1 indicate that /3, while 
clearly less than 2.5, differs significantly from for the cosmological models discussed in §4. The 
best parametric fit for the Einstein - de Sitter cosmological model (f3 = 0.84) is shown in Figure 1 
along with the expected relation for standard candle sources (j3 = 2.5). 

We now turn to a determination of the evolution of the luminosity. 

5.2. The Luminosity Evolution g(z) 

We examine the two commonly used parametric forms for the luminosity evolution, the ex- 
ponential and power law forms. These two different forms emphasize the evolution in different 
regions of the luminosity — redshift {L — z) plane. A correct parameterization will have the same 
value for its parameters when applied to samples with different limits (i.e. different coverage of 
the L — z plane). This fact can be used to test a given parametric form. Table 2 summarizes the 
results described in the subsequent sections. Figure 2 gives an example of the variation of the test 
statistic r as a function of the evolution parameter k. The optimal value of k, i.e. the value which 
indicates that L Q and z are not correlated, is given by the condition r = and the 1 a range of 
this parameter is given by the condition |r| < 1. These values are shown in Figure 2. 

5.2.1. Evolution of the form gk{z) oc e kt ^ 

Figures 3a and 3b display the best values and the 1 a ranges for the evolution parameter k as- 
suming both £l\ = and fij. = cosmological models for the different samples. The Durham/ AAT 
data exhibit much stronger evolution than the other data, indicating that this form of evolution 
does not adequately describe the data in the entire L — z plane. The values of k found for most 
of the data sets are considerably less than the previously determined value k pa 7.5 by Caditz and 
Petrosian (1990), which was dominated by the Durham/AAT data. 



5.2.2. Evolution of the form gk'(z) oc (1 + z) 
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Sample 


N 


T\ 


Pi 


T2 


P2 


Pi 


02 


Durham/ A AT 


419 


-0.65 


48.71 


-0.45 


35.02 


-0.22 


-0.18 


LBQS 


871 


4.65 


99.99 


3.75 


99.98 


1.11 


0.92 


HBQS 


254 


1.37 


82.92 


1.27 


79.42 


0.71 


0.59 


LBQS and HBQS 


1125 


4.80 


99.99 


3.98 


99.99 


1.01 


0.83 


Combined Data 


1552 


3.98 


99.99 


3.63 


99.97 


0.84 


0.70 



Table 1: Correlation data for the Quasar Hubble Diagram. N is the number of data points in each 
of the data sets of §2. The correlation r and the probability value P for rejection of independence 
between B magnitude and redshift z are given. The first set of values (tl, Pi) are found using the 
one-sided method of Efron and Petrosian (1992) and the second set of values (t2, P2) are found 
using the general method for doubly truncated data presented in §3. The parameters (3± and /?2 are 
found with this general method by fitting to the equation B(z) = B — j3 \og{d\{z)K{z)) + constant 
for the two cosmological models Dm = 1, £l\ = and Om = 0.15, £l\ = 0.85 respectively. 



Sample 


TV 


k 


kmin 


kmax 


k' 


k' 


y 

""max 


Durham/AAT 


419 


8.72 


6.66 


10.07 


3.53 


2.57 


5.05 


HBQS 


254 


5.39 


— 00 


6.49 


3.20 


— 00 


3.94 


LBQS 


871 


4.28 


2.66 


5.17 


2.02 


1.24 


2.53 


Combined Data 


1552 


5.15 


4.36 


5.70 


2.58 


2.14 


2.91 



Table 2: Values of the evolution parameters k and k' for the different data sets. These parameters 
refer to luminosity evolutions gk{z) oc e kt ^ and gk'(z) oc (1 + z) k \ respectively. The minimum and 
maximum values are those allowed at the 1 a level. Note the variation in k and the near constancy 
of k'. 
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g k -(z) = (i+z) k , n k = n A = o 

6 
3 



t- 



-3 
-6 

2 4 6 

k' 

Fig. 2. — A determination of the luminosity evolution parameter for the parametric form gk'(z) oc 
(1 + z) k ' . The correlation statistic r is shown as a function of k' for the combined data set for the 
Einstein - de Sitter cosmological model. We use the correlation test for multiply truncated data 
of §3, normalized so that r has a standard deviation of 1. The solid line at r = demonstrates 
the determination of the optimal value k' = 2.58 and the dashed lines at |r| = 1 demonstrate the 
determination of the 1 a region [2.14,2.91]. 
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o) g k (z) = e kt(z) , n A = o 

101 ' ' ' * ' ; ' i 
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_s □ 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 



b) g k (z) = e kt(z) , n k = 







A A A 



_X Combined Data 
Durham/AAT 
A LBQS 
a HBQS 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 



Fig. 3. — Variation of the evolution parameter k with fi^ for evolution of the form gk(z) cx e ht ^ 
for the three major surveys and for the combined data set. One a error bars are plotted for the 
CIm = 1 case. Error bars are similar for the other values of fi^- Note the significant difference 
between the deeper Durham/AAT data and the LBQS and HBQS data. This may indicate that 
the choice of gu{z) is not correct. Figure a) shows results for cosmological models with f^A = and 
Ofc = 1 — Qm- Figure b) shows results for cosmological models with 0,^ = and 17a = 1 — Qm- 

Figures 4a and 4b give results for the evolution parameter k' for the same cosmological models 
as above. The allowed ranges obtained from the different samples are much closer, thus this form 
of evolution is shown to be a closer approximation to the actual evolution than the exponential 
form. For the Einstein - de Sitter model, the best value of evolution parameter is k 1 = 2.58 with a 
one a range of k' 6 [2.14, 2.91]. As with the previous case, the values of k' are somewhat less than 
the previously found value of k' « 3 of Caditz and Petrosian (1990). 

It is clear that a better fit can be achieved with a different functional form for g{z) with two or 
more parameters. However, in what follows we assume the simpler form of evolution g(z) cx {l + z) k ' 
with k'{Q) the optimal value of k' for a given cosmological model Q. We therefore transform the 
data to {[L Q i, Zi]}f =1 with L (z,m,Cl) = L(z,m,Cl)/(l + z) k ' and apply the method of §3 to find 
non-parametric estimates for the density functions p{z) and ip(L ). This method now gives directly 
the cumulative functions a(z) = J 2 p(z)^dz and $>{L ) = f£? ip(L')dL' . 



5.3. The Density Evolution p(z) 



The cumulative density function <r{z) is the total number of objects within the angular area of 
the survey up to redshift z. If there is no density evolution, i.e. p{z) = p is a constant, then a cx V 
where V(z) is the co-moving volume up to redshift z. We determine a(z) and p(z) using the new 
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a) g,(z) = (1+z) k ', Q A = 



5 
4 

3 

2 

1 


0.0 0.2 0.4 0.6 0.8 1.0 1.2 




_^ Combined Data 
<> Durham/AAT 
A LBQS 
D HBQS 



b) g,(z) = (1+z) k ', n k = 




_X Combined Data 
Durham/AAT 
A LBQS 
a HBQS 



F ... i ... i ... i ... i ... i 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 



Fig. 4. — Same as Figure 3 with evolution of the form g^{z) oc (1 + z) k ' . The error bars nearly 
overlap at the optimal value for the combined data, indicating that this is an better parameterization 
over the observed range of the L — z plane. 

method for doubly truncated data. In order to determine if density evolution exists we fit a to V 
by a simple power law a(z) oc V x , where A / 1 indicates the presence of density evolution. If the 
density increases with redshift we expect A > 1 and if the density decreases with redshift we expect 
A < 1. Even if A = 1, however, density evolution may be present: the density may increase and 
decrease in such a way as to cancel and give a fit of A = 1. Figures 5a and 5b show the variation of 
a and p with V for the combined sample for three different cosmological models. The dotted lines 
show the best fits to the form a cx [V(z, For the Einstein - de Sitter model (top curves in 

Figure 5) we have A = 1.19, indicating that the co- moving density increases with redshift roughly 
as a x V 119 . This would indicate a simple power law density evolution p oc V 0Ag . The density 
evolution shown in Figure 5b exhibits this average behavior, but in detail is more complex: the 
density increases more rapidly at low z, reaches a plateau at z ~ 2 and possibly decreases at higher 
z. This behavior will be discussed below in more detail and for higher redshift data. 

As mentioned in §1, one cannot determine the evolution of sources (the function p{z)) and 
the evolution of the universe (the parameters f2j) simultaneously. Given one of these (e.g. the 
cosmological model), the other (the density evolution) may be determined from the data. The 
variation of A = dlncr/dln V with Um for the two different classes of cosmological models is shown 
in Figures 6a and 6b. It is evident that there is a monotonic variation of A with Om- In particular, 
A = 1 for the two cosmological models: S1m ~ 0, 12 a ~ 1 and ~ 0.15, S7a ~ 0.85. This second 
set of parameters is quite close to those currently favored by many observations. As above, this 
does not imply the complete absence of density evolution. The lower two curves in Figure 5b show 
the variation of p for these two models. Clearly, there is less variation than for the Einstein - de 
Sitter model, but the general behavior is similar. 



-18- 




8 10 12 9 10 11 12 13 

log(V - V min ) log(V - V min ) 



Fig. 5. — A plot of a) a = Jq p(z')dz' and b) p as a function of co-moving volume V(z, f2) for the 
combined data set for three different cosmological models. Curves i), ii) and hi) show results for 
Einstein - de Sitter, = 0.15, S7a = 0.85 and Qm = 0, f^A = models, respectively. Vertical 
normalizations are arbitrary. The straight dotted lines show the best fits to the parametric models 
a cc V x , with A = 1.19 for the Einstein - de Sitter model and A = 1 for the other two models. 
It is clear from b), however, that the model a oc V x is only a rough approximation. In fact, p 
undergoes more complicated evolution, increasing and peaking before decreasing at higher redshift. 
The luminosity evolution gk'(z) oc (l + z) k ' , with the optimal value of k' = 2.58, was used to remove 
the correlation from the original data. 
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Fig. 6. — The best fit for the power law index A = dlna/dln V is shown as a function of Qm for two 
different types of cosmological models. The index A is defined by the equation <j{z) oc V A , where 
V(z, CI) is the co- moving volume of space contained within the sphere of radius R(z). A value of 
A < 1 corresponds to a decrease in density with redshift and A > 1 corresponds to an increase in 
density with redshift. Pure luminosity evolution, A = 1, (a oc V(z)) occurs for models without 
cosmological constant at Om ~ ~ and for models without curvature at Qm ~ 0.15, ft a ~ 0.85. 
Figure a) shows results for cosmological models with VIa = 0. Figure b) shows results for models 
with = 0. 
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To further analyze the variation of p(z), in Figure 7 we show our non-parametric determination 
of p(z) for the Einstein - de Sitter cosmological model. Two sets of results are given. The first of 
these (depicted by squares) shows p(z) in the region 0.3 < z < 2.2 from the combined data. The 
second (triangles) shows p{z) in the region 0.3 < z < 3.3 from the LBQS data (which do not have 
a high redshift cutoff at z = 2.2) alone. As evident in Figure 7, the density increases relatively 
slowly at low redshift (p ~ (1 + z) 2 - 5 ) before reaching a peak at z « 2 and decreasing rapidly 
(p ~ (1 + z)~ 5 ) at higher redshift. As discussed further in §5.5, the decrease in density present in 
this data at redshift of about 2 is in agreement with high redshift survey results (Schmidt et al. 
1995, Warren et al. 1994). 



g,(z) = (1+z) k ', fi k = Q A = C 
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Fig. 7. — The non-parametric results for the density function p{z) for the Einstein - de Sitter 
cosmological model. Clearly, a simple model p oc (1 + z)^ will not adequately describe the density 
evolution. The squares show the results for the combined data in the region 0.3 < z < 2.2 and the 
triangles show the results for the LBQS data in the region 0.3 < z < 3.3. The density peaks at 
z « 2. 



5.4. The Luminosity Function ip{L ) 



In a similar fashion we may obtain the cumulative luminosity function &{L ) from the uncor- 
rected data set {L Q ,z}. Figures 8a and 8b show <J>(L ) for the combined data set with k' = 2.58, 
along with the best fits to a double power law form 

HLo) = (L /L.)k*+(L /L.)>»- (28) 
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We used the cosmological models Qm = 1, = (Einstein - de Sitter model) and Qm = 0.15, 
r^A = 0.85 (pure luminosity evolution with cosmological constant). In both cases, the results for the 
combined samples exhibit roughly double power law dependence on L with similar values of the 
fitting parameters. The primary differences between the two models are that the Einstein - de Sitter 
model gives a slightly gentler slope above the break luminosity, = 3.17 as opposed to = 3.59, 
and a lower break luminosity, = 6.19 x 10 29 erg / (sec Hz) as opposed to L* = 9.48 x 10 29 erg 
/ (sec Hz). 

a) g,(z) = (l+z) k ', n M = 1, n k = b) g,(z) = (l+z) k , fi M = 0.15, fi A = 0.85 




log(L ) log(L ) 

Fig. 8. — The cumulative luminosity function &(L Q ) = tp(L')dL' for the combined data set is 
shown as a function of L a = L/(l + z) k for two different cosmological models. The dotted curve is 
the best fit to a double power law model, equation (28). Figure a) shows results for the Einstein 
- de Sitter model. In this case, the best fit to the double power law form has break luminosity 
= 6.19 x 10 29 erg / (sec Hz) and power law indices k\ = 1.05 and ki = 3.17. Figure b) shows 
results for a cosmological model with £Im = 0.15 and f^A = 0.85. The best fit for this model has 
break luminosity = 9.48 x 10 29 erg / (sec Hz) and power law indices k\ = 1.16 and ki = 3.59. 

We check for possible variation in the shape of i^(L ) by dividing the data into three redshift 
bins: 0.3 < z < 0.86, 0.86 < z < 1.48, and 1.48 < z < 2.2. We then find the differential 
luminosity function ip(L ) for these three redshift bins, first assuming no luminosity evolution 
(g(z) = constant) and then assuming luminosity evolution g(z) oc (1 + z) k ' . These luminosity 
functions (with arbitrary vertical normalization) are shown in Figures 9a and 9b, respectively. The 
presence of a strong shift to higher luminosities is clearly evident for g(z) = constant. However, 
when the evolution g{z) oc (1 + z) k is taken out the luminosity function seems to exhibit little 
variation; the slopes appear roughly the same at low L Q and high L Q and the break luminosity does 
not vary as much with redshift. Although imprecise, these results indicate that our choice of g{z) 
removes most of the variation of the parameters «j with redshift, i.e. the shape of the luminosity 
function is almost invariant. 
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Fig. 9. — The luminosity function ip(L ) for three different redshift bins (containing equal numbers 
of sources) for the Einstein - de Sitter cosmological model. Vertical normalization is arbitrary. 
Figure a) shows results assuming no luminosity evolution, g(z) = constant. Figure b) shows results 
assuming luminosity evolution g(z) oc (1 + z) k ' with the optimal value of the parameter k' = 2.58 
found above. The second figure shows small variation in the shape of tp(L ) with redshift, indicating 
that most of the variation of the parameters ai is removed by the choice of luminosity evolution 
g{z) oc (1 + z) k ' . As expected, the high luminosity end can be fitted to a power law with index 
k 2 + 1 ~ 4. 
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5.5. The Luminosity Density C(z) 

Finally, we determine the luminosity density function C{z). This quantity is defined as the 
total rate of energy production by quasars in the optical range as a function of redshift; C{z) = 
Lip(L)dL. If the shape of the luminosity function is invariant then C(z) oc p(z)g(z). This rate 
depends in a complicated way on the distribution of masses of the central black holes and the 
variation of the accretion rate, both of which are related to the formation of galaxies and their 
evolution through mergers or collisions. Using the above results we can evaluate C{z) up to a 
redshift of 2.2. We extend this further to a redshift of 3.3 using the LBQS data, which is claimed 
to be complete up to this redshift. These results, with arbitrary vertical normalization, are shown 
in Figure 10. We first note the good agreement in the z < 2 region, indicating that perhaps the 
LBQS result at higher redshift is a representative behavior. We may also use high redshift surveys 
of quasars (Schmidt et al. 1995, Warren et al. 1994) to study C(z) in this range. Unfortunately, 
the selection of high z quasars in these samples is more complicated and the subsequent analyses 
involve more assumptions. For example, Schmidt et al. use the V/V max method to determine p(z), 
tacitly assuming that g(z) = constant (as well as ai = constant) so that C{z) oc p(z). We show 
these results (again with arbitrary vertical normalization) in Figure 10. These results agree with 
the general trend of decline in C(z) at high redshifts. 

It has been claimed (Cavaliere and Vittorini 1998, Shaver et al. 1998) that this rise and fall of 
C{z) with redshift is similar to the behavior of the star formation rate (SFR), which has recently 
been extended to high redshift (see, e.g. Madau 1997). We have shown this rate in Figure 10 as well. 
Although the general trend of rise and fall of the SFR and C{z) is the same, there is considerable 
difference in the detailed variation. The similarity may indicate some relation between the SFR 
and the feeding of the central engine of the quasars (e.g. both are affected by mergers). However, 
considering the many differences between star formation and the generation of energy by quasars, 
the observed difference between the SFR and C{z) in Figure 10 is not surprising. 

6. SUMMARY AND CONCLUSIONS 

Although there have been several analyses of quasar evolution in the past, our results differ 
from these previous results in two important respects: 

• We have used non-parametric statistical methods for multiply truncated data that allow us 
to combine samples with different selection criteria. 

• We have used the data to study models of the luminosity function that take into account 
both luminosity evolution g(z) and density evolution p{z). 

The new non-parametric statistical methods differ from those used in the past in the following 

ways: 




Fig. 10. — The luminosity density C{z) is shown as a function of redshift for the Einstein - de 
Sitter cosmological model. The squares and triangles give our results for the combined data and 
the LBQS data, respectively. The high redshift results of Schmidt et al. (1995) and Warren et al. 
(1994) are given by the letters "s" and "w", respectively. The vertical normalizations are arbitrary. 
All of the above results indicate that C peaks somewhere in the region z ~ 2. The stellar formation 
rate (SFR) as a function of redshift, as found by Madau (1997), is given by the asterisks. Although 
the SFR also exhibits the same general behavior as C(z), the peak occurs at a lower redshift and 
appears somewhat broader. 
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• No binning is required and most of the characteristics of the distribution functions are 
determined non-parametrically. 

• The methods are not limited to simply truncated data such as that found in flux limited 
surveys and can account for selection biases in generally truncated data where each data point has 
different truncation limits. In particular, these methods can treat samples with both upper and 
lower flux limits and redshift limits. 

• This versatility allows one to combine data from different surveys with different selection 
criteria. 

• The first of our techniques, a generalized non-parametric test of independence, allows one 
to determine the degree of correlation between luminosity and redshift, giving an indication of the 
luminosity evolution in the luminosity function. The evolution may then be determined paramet- 
rically. 

• The second of our techniques provides a non-parametric estimate for the univariate distri- 
butions in redshift and luminosity, i.e. the co-moving density evolution and the local luminosity 
function, respectively. 

We have applied these methods to the combined data from several large surveys and determined 
the luminosity evolution g(z), the density evolution p(z) and the luminosity function ip{L a = 
L/g(z)) of the generalized luminosity function of equation (^) for flat (Qk = and = 1 — Q\) 
and zero cosmological constant = and £Im = 1 — ^k) cosmological models. We assume 
a shape invariant luminosity function, on = constant. More complex luminosity functions, ^ 
constant, or those with luminosity dependent density evolution, etc., can be tested if the simpler 
prescription used here is not consistent with all of the data. We found that the scenario of equation 
@ provides an adequate description of the existing data. 

Our results may be summarized as follows: 

• We found a strong correlation between luminosity and redshift, indicating the presence of 
rapid luminosity evolution. 

• The parametric model of luminosity evolution (1 + z) k ' provides a better description of the 
data than the model e kt<yZ \ although neither parameterization perfectly removes the correlation 
in all areas of the L — z plane. In order to better model this evolution future analyses of quasar 
evolution could consider other parametric forms, including those with more than one free parameter. 

• The cumulative co- moving density of quasars may be approximately modeled as c(z) oc V x , 
where the value of A depends on the cosmological model. For example, A = 1.19 for the Einstein 
- de Sitter model and A = 1 for the cosmological models with VLm ~ 0, £lk ~ 1 and Om ~ 0.15, 
JIa ~ 0.85. This simple parameterization a oc does not describe precisely the variation of p with 
redshift. When examined in greater detail, the co-moving density shows a relatively slow increase 
(p ~ (1 + z) 2 5 ) for low redshifts and a rapid decline (p ~ (1 + z)~ 5 ) for z > 2 (for the Einstein 
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- de Sitter model). This is in qualitative agreement with the observed density evolution of high z 
quasars (Schmidt et al. 1995, Warren et al. 1994). Qualitatively similar behavior is found even 
for models that show no overall density evolution (A = 1). A more rigorous comparison cannot be 
made at this stage because the analysis of the high z data ignores possible L — z correlation and 
assumes pure density evolution. 

• The cumulative local luminosity function <&(L ) = J™ ip{x)dx has a double power law form. 
In the Einstein - de Sitter model the break luminosity is L* = 6 x 10 29 erg / (sec Hz) and the low 
and high luminosity power law indices are k\ = 1.05 and k<i = 3.17. There appears to be little 
variation with redshift of the shape of the cumulative and differential luminosity functions, thus 
the cti = constant prescription seems adequate. With more data one could determine precisely the 
variation with redshift of the shape of (fi(L ). 

• The above description of the luminosity function allows us to determine the rate of energy 
generation per unit co-moving volume of quasars as a function of redshift. We show that this 
function C(z) oc p{z)g{z) increases rapidly with z at low redshift, peaks around z ~ 2 and then 
decreases. This is also in rough agreement with the high z survey results mentioned above. This 
variation of C{z) is similar to but significantly different from recent determinations of the star 
formation rate. 

We would like to thank Bradley Efron for his invaluable help with the statistical methods used 
in this paper, Paul Hewett for help with the analysis of the LBQS data and Nicole Lloyd for helpful 
discussions. A. M. would like to acknowledge support from the Stanford University Undergraduate 
Research Opportunity program. 
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